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Abstract We combine a convectively driven dynamo in a spherical shell with a 
nearly isothermal density-stratified cooling layer that mimics some aspects of a 
stellar corona to study the emergence and ejections of magnetic field structures. 
This approach is an extension of earlier models, where forced turbulence simula- 
tions were employed to generate magnetic fields. A spherical wedge is used which 
consists of a convection zone and an extended coronal region to « 1.5 times the 
radius of the sphere. The wedge contains a quarter of the azimuthal extent of the 
sphere and 150° in latitude. The magnetic field is self-consistently generated by 
the turbulent motions due to convection beneath the surface. Magnetic fields are 
found to emerge at the surface and are ejected to the coronal part of the domain. 
These ejections occur at irregular intervals and are weaker than in earlier work. 
We tentatively associate these events with coronal mass ejections on the Sun, 
even though our model of the solar atmosphere is rather simplistic. 

Keywords: Magnetic fields, Corona; Coronal Mass Ejections, Theory; Interior, 
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1. Introduction 

Rece nt observations of the Solar Dynamics Observatory (SDO: [Pesnell. Thompson, and Chamber linl .[ 



2012| ) have provided us with a record of impressive solar eruptions. These erup- 
tions are mostly associated with coronal mass ejections (CMEs). These are events 
through which the Sun sheds hot plasma and magnetic fields from the corona into 
the interplanetary space. The energy causing such huge eruptions is stored in the 
magnetic field and can be released via reconnection of field lines (|Sturrock , Il980l: 
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Antiochos. De Vore. and Klimchukl . Il999l) . Some of the CMEs are directed to- 
wards the Earth, hitting its magnetosphere and causing phenomena hke aurorae. 
Furthermore, encounters with CMEs can cause sudden outages of GPS signals 
due to ionospheric scintiUation. The resuhing radiation dose from such events 
poses risks to astronauts. This is now also of concern to airlines, because the 
radiation load during polar flights can reach annual limits, especially for preg- 
nant women. This leads to great interest of scientists in many fields of physics. 
However, there is an additional motivation which comes along with space weather 
effects. The solar dynamo, which is broadly believed to be responsible for the gen- 
eration of the solar magnetic fie ld, needs to be sustained by s heddi ng magnetic 
helicity from the Sun's interior ( Blackman and BrandenbureL 2003 ). Mean- field 
and direct numerical simulations have shown that the magnetic field generation is 
catastrophically quenched at hig h magnetic Reynolds numbers in closed systems 
( Vainshtein and Cattaneol 19921) that do no t allow magnetic helicity fluxes ou t of 
the domain ( Blackman and Field. l2000allbl: [Brandenburg and Sandinl.[200l. or 
between different parts of it |Bra ndcnburg. Cand elaresi. and Chatteri ^. 120091: 
Mitra et al\ . \20ld iHubbard and B ra ndenbur g, The magnetic Reynolds 



number, which quantifies the relative importance of advective to diffusive terms 
in the induction equation, is known to be very large in the Sun, therefore 
implying the possibility of catastrophic quenching in models of the solar dy- 
namo, unless efficient magnet i c heli city fluxes occur, for example through CMEs 

( Blackman and Brandenburg . 20031). Indeed , CMEs are well known to b e closely 

associated with magnetic helicity (jLowl, 2001 ) . In particular obs ervations ( Plunkett et al\M 



20001: 



Regnier. Amari. and Kersalel 20021) and a recent study bv lThompson. Kliem. and Torok^ 



(j201l[ ). where the observations are compared with numerical models, suggest 
that CMEs have a twisted magnetic structure, implying that CMEs transport 
helicity outwards. 

There has been significant progress in the study of CMEs in recent years. 
In addition to improved observations from s pacecrafts, e.g . SDO or the Solar 
TErestical RElation Observatory f STEREO: iKaiser et ai\ . l2008l )). there have 
also been ma j or advances in the field of n umerical modeling of CME events 
( Russev et al. . l2003t lArchontis et al. I, Eooi). However, the formation and the 
origin of eruptive events like CMEs is not yet completely understood. Simulating 
CMEs and their formation is challenging. Leaving the difficulties of modeling the 
interplanetary space aside, a CME, after being ejected into the chromosphere or 
the lower corona, travels over an extended radial distance to the upper corona. In 
this environment, density and temperature vary by several orders of magnitude, 
which is not easy to handle in numerical models. Additionally, the origin of the 
CMEs is assumed to relate to the magnetic fields and the velocity pattern at 
the surface. However, the surface magnetic and velocity fields are rooted in the 
solar convection zone, where convective motions, in interplay with differential 
rotation, generate the magnetic field and the velocity patterns that are observed 
at the surface. The majority of researchers modeling CMEs do not include the 
convection zone in their setup, and thereby neglect the effect of the magnetic 
and velocity fields being rooted to this layer. Often the initial conditions for the 
magnetic and velocity fields are prescribed or taken from 2D observations; see 
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for example Antiochos. De Vore. and Klimchuk ( 1999t ) and Amari et al. ( 1999t ) 
as well as Torok and Kliem ( 2003[) . 

Another approach is to study the emergence of flux ropes from the lower con- 
vection zone into the corona. In the presence of stro ng shear, convection s i mula- 



tions have been sho wing the formation of flux tubes (jGuerrero and Kapvlal . 12011 



Nelson et al\ . 120111 ) , but such structures are similar to vortex tubes whose diame- 



ter is known to relate with the visco-resistive scale (jBrandenburg. Procaccia. and SegelL l 



19951 ) ■ In other approaches flux ropes are inserted in a self-consistent model, but 



their origin is left unexplained. In several r ecent papers ( Martmez-Svkora. Hansteen. and Carlssonl .[ 
20081: IJouve and Brunl . [2009t iFang et the focus lies on the emergence 



of magnetic flux and the resulting features in the solar atmosphere. However, 
eruptive events have not been inv estigated with this setup. In earlier work 
( Warnecke and Brandenburg . 2010l . hereafter WB) a different approach was de- 
veloped. The solar convection zone was combined with a simple model of the solar 
corona. The magnetic field, which was here generated by dynamo action beneath 
the solar surface, emerged through the surface and was ejected out of the domain. 
The focus was on the connection of the dynamo-generated field and eruptive 
events like CMEs through the dynamo-generated twist. WB used a simplified 
coronal model and drove the dynamo with forced turbulence. These simplifica- 
tions allowed them to study the emergence and a new mechanism to drive ejec- 
tions in great detail. In subsequent work (jWarnecke. Brandenburg, and Mitral . 



20111 . hereafter WBM), the setup of WB was improved by using a spherical 



coordinate system and helical forcing with opposite signs in each hemisphere to 
mimic the effects of rotation on inhomogeneous turbulence. In addition, WBM 
included the stratification resulting from radial gravity for an isothermal fluid. 
To improve this model, we now employ c onvec tion to generate the velocity field. 
In a related approach, iPinto and Brunl (|201l[ ) considered convective overshoot 
into the chromosphere and the excitation of gravity waves therein, but dynamo- 
generated twist seemed to be unimportant in their work. The turbulent motions 
driving the generation of magnetic field arc now self-consistently generated by 
convective cells operating beneath the surface. The setup o f the convection zone 
follow s id eas of Kapvla. Korpi. and Brandenburg ( 2008 ). iKapvla et'al. (|2010l . 
2011 ) and Kapvla. Mantere. and Brandenburg! ( 20 111 " 2012 ) . There are other ap- 



proaches simulati ng convection in hot m assive stars, which have thin subsurface 



convection zones (jCantiello et a/.l . l201ll ). But we now use an extended cooling 



layer to describe some properties of a solar corona. The results of this work 
complement those of earlier work and can be compared with observations. The 
model of the solar atmosphere is still a very simplified one, but can be regarded 
as a preliminary step, which will provide a reference point for improved work in 
that direction. 



2. The model 



As in WB and WBM, a two-layer model is used, which represents the convec- 
tion zone and an extended corona-lik e layer in one and the sam e model. Our 
convection zone is similar to those of Kapvla et al. ( 2010l 2011 ). The domain 
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is a segment of the Sun and is described in spherical polar coordinates (r, 0, (p). 
We model the convection zone starting at radius r = 0.7 i? and the solar corona 
until r = Rc, where Rc = 1.5 R in the present models, where R corresponds to 
the solar radius. In the latitudinal direction, our domain extends in colatitude 
from 6 = 15° to 165° and in the azimuthal direction from (p ^ 0° to 90°. We 
solve the following equations of compressible magnctohydrodynamics: 



dA 

'dt 
Dlnp 



U X B + T]V'A, (1) 



Dt 

— = g-2fnoxU+-(JxB-Vp + V- 2i/pS) - D{r, 9, t), (3) 
Dt p 

= -V.A'VT + 2:.S2 + ^j2-reooi, (4) 
Dt p p 

where the magnetic field is given by S = V x A and thus obeys V • B = at 
all times, /xq is the vacuum permeability, 77 and v are the magnetic diffusivity 
and kinematic viscosity, respectively, D/Dt = d/dt + U ■'V is the advective time 
derivative, p is the density, and U is the velocity. The traceless rate-of-strain 
tensor is given by 

where semicolons denote covariant differentiation; see Mitra et al. ( 20091 ) for 



details. = fig (cos 0, — sin 0,0) is the rotation vector, p is the pressure, K is 
the radiative heat conductivity, and D{r,6,t) describes damping in the coronal 
region; see Section for details. The gravitational acceleration is given by 

g = -GMr/r^, (6) 

where G is Newton's gravitational constant, and M is the mass of the star. The 
fluid obeys the ideal gas law, p = (7 — l)pe, where 7 = Cp/cy = 5/3 is the 
ratio of specific heats at constant pressure and constant volume, respectively, 
and e = CyT is the internal energy density, which defines the temperature T. 
The cooling term Fcooi will be explained in Equation (jlO|) below in more detail. 

2.1. Initial setup and boundary conditions 

For the thermal stratification in the convection zone, we consider a simple ana- 



lytical setup instead of profiles from solar structure models as in, e.g.. lBrun et al 



(|2004[ ). The hydrodynamic temperature gradient is given by 

dT -\q\ 



dr c„(7 - l)(m + 1) ' 



(7) 



where m = m{r) is the radially varying polytropic index, for which we assume 
a stepwise constant profile. We also use Equation (O as the lower boundary 
condition for the temperature. This gives the logarithmic temperature gradient 
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Figure 1. Initial stratification of temperature (dashed line), density (solid), pressure (dot- 
dashed) and the Brunt- Vaisala frequency A^^ = —{\g\/ Hp){^ — V^d) (dash-triple-dotted) for 
Run A5. The subscripts h refers to the values at r = 0.7 i?. The dotted horizontal (vertical) 
line denotes the value of zero {r = R). 



V (familiar to those working in stellar physies, but not to be confused with the 
operator V) as 

_ d\nT _ 1 

^ d\np ^ m + l ^ ' 

The stratification is convectivcly unstable if V — Vad > 0, where Vad = 1^1/7 
is the adiabatic temperature gradient, corresponding to m < 1.5 for unstable 
stratification. We choose m = 1 in the convectively unstable layer beneath the 
surface, r < R. The region above r = i? is stably stratified and isothermal 
due to a cooling term Fcooi with respect to a constant reference temperature 
in the entropy equation. The density stratification is obtained by requiring the 
hydrostatic equilibrium condition to be satisfied. 

The thermal conductivity follows from the constancy of the radial luminosity 
profile L{r) = Lq = const throughout the domain and is given by 

K ^ (9) 

To speed up the thermal relaxation processes, we apply shallower profiles, corre- 
sponding to p oc T^ *, for the thermal variables within the convectively unstable 
layer. The value m = 1 is just used in the convection zone to determine the ther- 
mal conductivity. In Figure [1] we show the initial non-convccting stratification. 
The radial temperature gradient at the bottom of the domain is set to a constant 
value, which leads to a constant heat flux into the domain. In the coronal part 
the gradient goes smoothly to by using the r dependent cooling function Fcooi, 
which is included in the entropy evolution The cooling term is given by 

Fcooi = Fo/(r) (^^) , (10) 



SOLA: sol_paper.tex; 2 October 2012; 4:58; p. 5 



J. Warnecke et al. 



where /(r) is a profile function equal to unity in r > R and smoothly connecting 
to zero in r < _R, and Fq is a cooling luminosity chosen so that the sound speed 
in the coronal part relaxes towards c^q = Cg(r = Rc). Whether the stratification 
is convectively stable or not depends on the Brunt- Vaisala frequency N, defined 
through 

where Hp — ~dr/d\np is the pressure scale height. If iV^ is negative, the 
stratification is unstable. 

We initialize the magnetic field as a weak, random, Gaussian-distributed seed 
field in the whole domain. In the coronal part the magnetic field diffuses after 
a short time. We do not use a background coronal field, so the field is self- 
consistcntly generated by the dynamo in the convective layer. We apply periodic 
boundary conditions in the azimuthal direction over a 90° fraction of the full 
circumference. For the velocity we take stress-free boundary conditions on all 
other boundaries. As in WBM, the stress-free boundary conditions prevent mass 
flux, so no stellar wind is possible. Because no mass can escape, material will 
eventually fall back from the boundary. Thermodynamic variables have zero 
gradients at the latitudinal boundaries. We employ perfect conductor boundaries 
for the magnetic held at the latitudinal and at the lower radial boundaries, and 
radial field conditions at the outer radial boundary. The latter is motivated by 
the fact that in the Sun, the solar wind pushes the magnetic field to open field 
lines and at a radius of r = 2.0 . . . 2.5 solar radii, the field lines are mo s tly ra - 



dial (jLevine. Schulz. and Frazierl . 119821 : iHoeksema. Wilcox, a nd Schcrrer', '^198 



This c hoice has been substantiated by subseq uent work of iWang and Shcele 



jelevl 



(jl992n as well as lSchriiver and De Rosal (j2003f ). While this choice might still be 



too restrictive for coronal holes and coronal streamers, and given also that our 
radial extent in most of the simulations is smaller than r = 2R, we nevertheless 
choose the vertical field boundary condition because it satisfies our primary 
objective of letting magnetic helicity leave the domain, w hich is believed to be 
crucial f or the dynamo to operate at lar ge values of Rej\,f (jBlackman and Fieldl . 



2000al lb: Brandenburg and Sandinl . 20041 ). However, we must be aware of the fact 



that with this choice our description of the field in the exterior layer is not a 
realistic one. 

To describe the corona as an isothermal extended cooling layer is a serious 
simplification, in that the temperature inside the coronal layer is not higher 
than in the convection zone as in a real stellar corona, but it stays fixed at the 
surface value; see Figure [1] Besides the fact that a simple cooling layer is easy to 
handle numerically, we emphasize the importance of facilitating comparison with 
previous models of WBM. It can also be seen as a step towards studying effects 
that are not solely due to a low plasma /3 corona, for which the magnetic pressure, 
i.e. the magnetic field, is strong compared with the gas pressure (/3 = 2^op/B^). 
Indeed, given that our initial field is weak, the plasma f3 is necessarily large in the 
outer parts. We note that it is not even clear whether a hot corona promotes or 
hinders coronal ejections. To understand the formation and evolution of magnetic 
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ejections, studies that isolate these effects, such as the present one, may be 
important. 

We use the Pencil Code0 with sixth-order centered finite differences in space 
and a third-order accurate Runge-Kutta scheme in time; see iMitra et al. I (l2009l) 
for the extension of the Pencil Code to spherical coordinates. We use a grid 
size of 128 x 128 x 64 mesh points (Runs A5 and Arl), and 256 x 256 x 128 
(Run A5a). 

2.2. Velocity damping in the corona 

Whether the solar corona rotates like a solid body or differentia lly coupled with 
the photosphere is unclear. In recent work bv lWohl et all (l20ld) . where SOHO- 
EIT data of the bright points in the solar corona were used to estimate the 
rotation speeds, it was found that the corona rotates similarly as the small 
magnetic featur es in the photosphere. Similar results have been obtained by 
Badalyan ( 2010l ). where the coronal rotation has been measured by analyzing 
the green Fc''^^^ 530.3 nm line. This author finds also a variation pattern with the 
activity cycle. However, the observations of the "boot" coronal hole by SKYLAB 
suggested rigid rotati on (jTimothv et all 119751 ) . Recent work on coronal holes by 



Lionello et al\ ( 2005 ) claims that the rigid rotation is only an apparent one. 
The magnetic field is sheared by the differential rotation, but the boundary of 
the hole remains relatively unchanged, due to rcconncction. Owing to the low 
plasma j3 in the solar corona, the fluid motions are dominated by the magnetic 
flelds whose footpoints are anchored in the photosphere or even further down. 
So the magnetic field might then be rigid enough to prevent differential rotation 
of the solar corona. However, the observed bright points and other features in 
the corona are strongly correlated with the magnetic field so they can give a 
misleading picture about the global rotation of the corona. 

In our simulations, the Coriolis force is included in the momentum equation 
as a consequence of the rotation. In the solar corona the density is more than 14 
orders of magnitude smaller than in the lower convection zone. Because of the 
weak density stratification in our simulation, the Coriolis force in our coronal 
part is too strong and can cause possible artifacts such as the magnetorotational 
instability. To avoid this — at least for runs with rapid rotation — we apply a 
damping function D{r, 9) in the momentum equation, which is given by 



D{r, 9, t) = — e(r - R) U{r, 9, t), 



where 



e(r - i?) = i 



1 + tanh 



r-R 



(12) 



(13) 



with Td being the damping time and w the width of the transition layer from 
convection zone to the coronal part. Here and elsewhere, the overbar denotes 
averaging over (f>, defined as F{r, 6, t)= J F{r, 9, (j), t) d^/27r. Occasionally we also 
use time averages denoted by (.)(. 



^http : //pencil-code . googlecode . com 
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Table 1. Summary of the runs. Rc is the fluid Reynolds number, Urms = y 3/2(?7^ + Ug) 

is the volume-averaged rms velocity in the convection zone normalized by ^/ GM/R, Ptm 
is the magnetic Prandtl number, Co is the Coriolis number, and h^^i is the maximum value 
of the relative kinetic helicity using azimuthal averages as defined in Section 12.31 ^ and 

^ give the density ratios of the bottom of the convection zone to those at the surface 
and the top of the domain, respectively. In the right-most column we note if damping for 



velocity 


in the coronal part is used (Y) or 


not (N)' 


; sec Section 12.21 










Run 


Resolution 


^rms 


Re 


PrM 


r>2 

rms 


Pb 


Pb 


Co 


hrcl 


D 


y/GM/R 




Ps 


Pt 


A5 


128^ X 64 


0.0072 


3.3 


10 


0.1-0.4 


3.6 


39 


7 


0.5 


N 


A5a 


256^ X 128 


0.0105 


100 


1 


0.2 


3.6 


39 


4.5 


0.3 


N 


Arl 


128^ X 64 


0.0040 


38 


1 


1.5-5.5 


3.6 


39 


50 


0.3 


Y 



2.3. Units, nondimensional quantities, and parameters 

Dimensionless quantities are obtained by setting 

R^GM = pb = cp = fio^l, (14) 

where pb is the density at r = 0.7 R. Below, we wiU describe the properties of 
the runs by the foUowing dimensionless parameters: fluid Reynolds number Re = 
u-n-as/vkf, magnetic Reynolds number Rca/ = Urms/?7^f, where kf = 27r/0.3i? is 
an estimate for the typical wavenumber of the energy-carrying eddies and Urms ~ 
■\/3/2([/2 _j_ [/2^ 13 lY^Q volume-averaged rms velocity in the convection zone (r < 
R). In our definition of Urms we omit the contribution from the (/(-component 
of the velocity, because it is dominated by contributions from the large-scale 
differential rotation that develops when rotation is included and would give an 
atypical estimate of the convective turnover time. To compensate for this, and 
to have an estimate of Urms comparable with earlier work, we apply the 3/2 
correction factor. We also define the magnetic Prandtl number Vim ^ i^/v = 
Rej\//Re and the Coriolis number Co = 2r2o/urmsfcf • Time is expressed in units 
of r = (urms^f) , which is the eddy turnover time in the convection zone. 
We measure the magnetic field strength as the rms value averaged over the 
convection zone i?rmsj where we often normalize this value with the equipartition 
value of the magnetic field defined by B^^ = /io(pM^jjjg)),.<fl. The relative kinetic 
helicity is h,-ci{r, t) = u) ■ w/wrmsUrmsj where u; = V x m is the vorticity and Wmis 
is its rms value inside the convection zone. 



3. Results 

3.1. Hydrodynamic phase of the simulations 

After around 100 turnover times, convection has reached saturation and we 
find convection cells as typical patterns in the radial velocity just below the 
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Figure 2. Flux balance from Run A5. The different contributions to the total luminosity 
(solid blue line) are due to radiative diffusion (dashed red line), resolved convection (dotted 
black line) and cooling (dash-triple-dottcd black line). The black thin line denotes the zero 
level and the surface [r = R), respectively. 



surface. In our model, the dominant energy transport mechanisms are radiative 
and convective fluxes in the bulk of the convection zone and an (optically thin) 
cooling flux in the outer (coronal) parts. The radiative and convective fluxes are 
defined as: 



-F'rad = -A'— , Fconv = Cppu'^T', (15) 

where the averages are taken over 9 and (f> and the prime indicates fluctuations 
about the respective mean quantity. In our present setup, however, the convective 
flux reaches barely about 5% in the convection zone; see Figure [2] where we plot 
the relevant contributions to the luminosity for Run A5. Above the surface the 
cooling takes over to maintain an approximately isothermal atmosphere. The 
total flux is constant, except for small departures near the surface. The kinetic 
energy and viscous fluxes are negligible in the present runs. 

To determine the degree of overshooting and penetration into the stably strat- 
ified layers above the convection zone, we show in Figure [3] the radial velocity 
above the surface at r = 1.15, 1.25, and 1.35 R for Run A5. At low latitudes, there 
is very little radial penetration (velocity features are only seen until r = 1.15 i?), 
while at higher latitudes the radial velocity pattern is transmitted all the way to 
1.35 i?. This is not surprising in view of the Taylor-Proudman theorem, which 
states that for rapid rotation (large values of Co) the local angular velocity of 
the gas is constant along cylindrical surfaces. 

Next, we plot in Figure 2] the rms values of all three velocity components for 
Run A5. The amplitude of the radial velocity component falls off the fastest. 
The latitudinal component also falls off with radius, but remains about three 
times larger than the radial component. The longitudinal component, on the 
other hand, increases with radius in a way that is compatible with rigid rotation 
with an angular velocity that is somewhat larger than the rotation rate of the 
frame of reference. 
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Figure 3. Radial velocity {Ur) above the surface for r = 1.15, 1.25, 1.35 R from left to right, 
for Run A5. Dark blue shades represent negative and hght yellow positive values. 



0.008 




0.000 t 



Figure 4. Root-mcan-square values of Ur (dotted), Ug (dashed), and l/^ (dash-dotted) as 
a function of radius for Run A5. The solid line shows the radial profile of our nominal rms 

velocity, Urms = \J'i/'i{Ur + ^g) ■ '^^^ (red) vertical line indicated the surface at [r = R). The 

values arc normalized by GM / R. 



The size of the convection ceUs depends st rongly on the strength of rotation 
and th e degree of density stratification; see also Kapvla. Mantere. and Brandenburg 
(|201l[ ). We plot the radial velocity Ur &t r = 0.89 R for Runs A5, A5a, and Arl 
in Figure [5] The Run A5 has a low fluid Reynolds number and therefore the 
convection cells are large; see Table [1] The flow pattern shows clear 'banana 



cells' as in previous work with comparable Coriolis parameter, cf. iKapvla et 



(|201l[ ). A higher fluid Reynolds number and higher resolution, as in Run A5a, 
allow the velocity field to form more complex structures. However, the banana 
cells are still visible. If one now looks at a simulation with more rapid rotation 
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Figure 5. Radial velocity (Ur) beneath the surface (r = 0.89 i?) for Runs A5, A5a, and Arl 
from left to right. Dark blue shades represent negative and light yellow positive values. 
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Figure 6. Differential rotation profiles f2(r, 6) = U ^/{r sin6) + Qq for Runs A5, A5a, and 
Arl from left to right. Dark blue shades represent low and light yellow high values, overplotted 
by the isocontours with solid black lines. The dotted white lines parallel to the rotation axis 
are given for orientation and the dashed line indicates the surface (r = R). 



(Run Arl, plotted in the right-most panel of Figure[5]) with a Coriolis number of 
Co = 50, the number of banana cells increases and they arc more clearly visible 
than in Run A5a. Note also that the radial velocity is now significantly reduced 
at high latitudes inside the inner tangent cylinder. 

In the Sun, differential rotation is an important element to produce the 
magnetic field structures observed at large scales, exhibiting a cyclic behavior 
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over time, as manifested by the sunspot cycle. To illustrate the differential 
rotation profiles generated in the simulations, we plot the azimuthally aver- 
aged angular velocity, fl{r,9) = U^/{rsm9) + fio, for Runs A5, A5a, and Arl 
in the saturated state of the simulation; sec Figure IHl In the plot, we show 
isocontours of angular velocity with solid black lines. In the convection zone 
the contours of local angular velocity tend to be cylindrical, which is likely a 

consequence of the absence of a strong latitudinal modulations of specific eiitropy 

j Brandenburg. Moss, and TuominenI Il992 ; Kitchatinov and Riidieer , 19951 : iMiesch. Brun. and Toomr3 .[ 

The coronal part seems to rotate as a solid body outside the outer tangent 
cylinder (i.e., for rsinO > R), while inside it some differential rotation occurs 
also in the coronal part. In the convection zone between the inner and outer 
tangent cylinders, the angular velocity is enhanced relative to that inside the 
inner tangent cylinder (sec the first and second panels of Figure 15]), while in the 
case of extremely rapid rotation this may actually be reversed. 

In the three runs shown in Figure [5] the stratification in the whole domain 
is just Ph/pt = 40, which is rather small compared to the stratification of 
the Sun {ph/pt ~ 10^*). It seems, therefore, that the Coriolis force is acting 
much more strongly in the coronal part of our simulation than in reality. In 
the Sun the Lorentz force plays a more important role in the corona than in 
our model. In the convection zone, we find quenching of convection due to 
rapid rotation. In Run A5, where Co = 7, the lines of constant rotation rate 
arc more radial than vertical and show super-rotation, i.e., the equator rotates 
faster than the poles. As expected, this tends to coincide with locations wher e 



Riidieer . 



the Reynolds stress in the radial direction is negative (see, e.g. 
However, the convection cells are rather big and have a strong local influence 
on Ucf, an d could in principle lead to sub r otatio n; see the corresponding dis- 



cussion m iDobler. Stix. and Brandenburd (|2006f ). Note that the rms velocity 
in Run A5 is two times smaller than in Run A5a, which has higher resolu- 
tion and higher fluid and magnetic Reynolds numbers (Rc ~ Re^/ ~ 100). 
Due to this, we find clear super-rotation, even though the Coriolis number 
is slightly lower (Co = 4.5) than what is realized in Run A5. In the third 
case. Run Arl, where the rotation is extremely rapid (Co — 50), we also find 
super-rotation, where the lines of constan t rotation rate are almost all vertical . 
In co mparable work ( Kapvla et al . l201lt iKapvla. Mantere. and Brandenbrirg . 



20111 ). super-rotation has been found, when the Coriolis number was larger 
than 4. This is similar to our results including a coronal part. In addition, 
there is a minimum of the rotation rate at mid-latitudes and a polar vortex at 
high latitudes. Rotation profiles, which show a comparable behavior, have been 
found by several groups (Miesch et al\. 2000t Elliot . Micsc h, and Toomre . I2OOOI: 



Kapvla et all 1201 ll : iKapvla. Mantere. and Brandenburg. 20111 ) . The region with 



the higher rotation rate near the equator is limited to the upper convection zone 
and can even penetrate into the coronal part. In Run Arl the velocity damping 
described in Section is used. By comparing the right-most panel of Figure [51 
with damping, to the left-most panels, without it, we conclude that the damping 
does not make much of a difference to the coronal velocity structures. 

Simulations with randomly forced turbulence (WB,WBM) have shown that 
the relative kinetic helicity hy^i has a strong influence both on the generation of 
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Figure 7. Relative helicity h,-^i{r,t) = lj ■ u/ajrms'irms plotted for Runs A5, A5a, and Arl 
from left to right. Dark blue shades represent negative and light yellow positive values. The 
dashed line indicates the surface (r = R). 



large-scale magnetic fields and the ejection events. In WB and WBM, values of 
hrei of order unity were achieved by using a forcing function with purely helical 
plane waves. In the convection runs presented here, however, values of large 
relative helicity, h-^d ~ 0.5, are obtained (for Run A5), at least at certain radii. 
In Figure [71 we present contour plots of azimuthally averaged relative helicity 
in the meridional plane for Runs A5, A5a, and Arl. All three show the typical 
sign rule of kinetic helicity under the influence of rotation, i.e. the northern 
hemisphere has predominantly a negative sign and the southern a positive one. 
Close to the bottom of the convection zone, the sign changes, which has earlier 



been reported by several authors both in Cartesian (e.g. [Brandenburg et al 
lOssendriiver. Stix. and Brandenburg . l200l[ ) and spherical geometries (e.g. 
Miesch et al. . I2OOOI: iKapvla et adlioiol) . Only in Run Arl with rapid rotation, 
the behavior is not that clear. The relative helicity is no longer confined to the 
convection zone, but significant values occur also in the coronal region. The sign 
rule still holds within the convection zone, while a more complicated sign behav- 
ior is visible in the coronal part. The maximal values of the azimuthally averaged 
helicity arc around /iroi = 0.3, occurring close to the surface. In Run A5a, the 
maximum value is slightly higher and is located in the middle of the convection 
zone, although relatively high values are present in the coronal part as well. 
It is not yet completely clear how high values of relative kinetic helicity can 
be achieved; strong rotation tends to suppress it, whereas strong stratification 
increases it. Its exact role in generating coronal ejections is yet unclear. 
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Figure 8. Time-averagod for Run A5, Ar5a, and Arl from left to right. Dark blue shades 
represent negative and light yellow positive values. The magnetic field is normalized by the 
equipartition value. The dashed line indicates the surface {r = R). 



3.2. Convective dynamo 

The convective motions generate a large-scale magnetic field due to dynamo 
action. The magnetic field grows first exponentially and begins then to affect the 
velocity field. The effects of this backreaction can be subtle in that we found a 6% 
enhancement of the rms velocity after saturation. The growth of the magnetic 
field saturates after around 200 to 1000 turnover times, depending on the Coriolis 
and Reynolds numbers. In the runs in Table [TJ we obtain different dynamo 
solutions for the saturated field. 

In Figure [5] we show the time-averaged azimuthal magnetic field for 
Run A5, A5a, and Arl. Note that the (j) component of the magnetic field is 
also strong in the coronal part and roughly antisymmetric about the equator. 
Furthermore we find an oscillation of the volume-averaged rms magnetic field in 
the convection zone; see the left-hand panel of FigurcElfor Run A5. The growth 
tends to be steeper than the decline, the period being around t/r = 220. The 
field reaches a maximum of 60% of the equipartition field strength, i?eq, which 
is comparable to the values obtained in the forced turbulence counterparts both 
in Cartesian and spherical coordinates (WB,WBM). Comparing this with the 
change of the kinetic energy, plotted as fiuctuations of the rms velocity squared, 
we find an anti-correlation with respect to the magnetic field oscillation. The 
magnetic field is high (low), when the velocity is low (high). In the work by 
Brun. Browning, and Toomrd (|2005l ). the authors interpret this behavior as the 
interplay of the magnetic backreaction and the dynamo effect of the differential 
rotation. Due to the Lorentz force a higher magnetic field strength leads to 
quenching of the differential rotation. An increased magnetic field quenches 
the Reynolds stress and thus lowers the differential rotation, which limits the 
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Figure 9. Phase relation of the magnetic field (B'^^^/ B^^, solid black lines) and the velocity 
field (urma(<)'^/(^'?ins)i> dashed red lines) in the convection zone for Runs A5 (left panel) and 
Arl (right). The velocity has been multiplied by a factor of 0.3 (left panel) and 3 (right), 
respectively, and smoothed over five neighboring data points. 




Figure 10. Variation of B^ and Br in the convection zone at r = 0.9R for Run A5. Dark blue 
shades represent negative and light yellow positive values. The dashed horizontal lines show 
the location of the equator at = 7r/2. The magnetic field is normalized by the equipartition 
value. 



magnetic field. A weak magnetic field leads to stronger differential rotation. 
Similar behavior has been observed also in previous forcing simulations (WBM). 
This behavior is not seen as clearly in the large-scale magnetic field which shows 
variations in strength, but not in sign. As shown in Figure fTUl for Run A5, the 
and Br have local maxima in time and in latitude, but the overall structure is 
nearly constant in time. Even though the large-scale field structure is stationary, 
the small-scale structures show an equatorward migration near the equator. The 
reason for this is unclear, but meridional circulation does not seem to play a role 
here. 
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Figure 11. Variation of and Br in the convection zone at r = 0.9R for Run Arl. Dark blue 
shades represent negative and light yellow positive values. The dashed horizontal lines show 
the location of the equator at 6 = n/2. The magnetic field is normalized by the equipartition 
value. 





Figure 12. Variation of in the convection zone at r = 0.9R for Run A5 (left panel) and 
Run Arl (right panel). Dark blue shades represent negative and light yellow positive values. 
The dashed horizontal lines show the location of the equator at 9 = it/2. The velocity is 
normalized by the mean rms velocity in the convection zone. 



In Run Arl, the magnetic field reaches up to 5.5 times the equipartition value, 
but does not show a periodic oscillation; see the right hand panel of Figure O 
In comparable work (|KaDvla et all |2010[ ). similar values for the field strength 
were found. However, the rms velocity is also quenched, when the magnetic field 
is high. Looking at i?^ and B^, plotted over time and latitude in Figure [TTl the 
large-scale magnetic field is similar to Run A5, which is constant in time with- 
out any oscillation. In the recent work by iKapvla. Mantere. and Brandenburg! 
(j2012n , the authors found an oscillatory behavior of and Br including equa- 
torward migration for latitudes below 60°, which is the first time that such 
a result is obtained from direct numerical convection simulation of rotating 
convection. 
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Figure 13. Time series of a coronal ejection near the equator {9 = 7r/2), taken from Run A5. 
The normahzed current helicity, fioR J ■ B / (B'^)t, is shown in a color-scale representation 
from different times; dark blue represents negative and light yellow positive values. The dashed 
horizontal lines show the location of the surface a,t r = R. 



The azimuthal velocity versus time and latitude fFigure fT^ shows minima 
at the same times as the maxima of the magnetic field occur. In Run A5a, the 
occurrence of strong magnetic fields suppresses the differential rotation. The 
pattern of the azimuthal velocity is symmetric about the equator and shows an 
oscillatory behavior, which is not that clear in the large-scale magnetic field. In 
the plot in Figure [T^] of Run Arl, we find just one localized minimum, which 
coincides with the low values of Urms(i)^/ (u?ms)t between t/r = 2100 and 2400. 

3.3. Coronal ejections 

In the runs that we have been performed so far, and of which only three have 
been discussed in this paper, only a small fraction of events can be identified with 
actual coronal ejections similar to the ones seen in WB and WBM. Especially 
the Runs A5 and Arl show some clear ejection events. There the magnetic field 
emerges out of the convection zone and is ejected as an isolated structure. In 
Figure [13] we have plotted the normalized current helicity, fioRJ-B/{B'^)t, as 
a time series for Run A5. At small scales, the current helicity density, J ■ B, 
is a good proxy for magnetic helicity density, A ■ B, and is, as opposed to the 
latter, gauge invariant. In addition, the current helicity can be an indicator of 
helical m agnet i c structures, which are bel i eved t o be present in coronal mass 



ejections (ILow . 1994 . 200ll Plunkett et al. . 2000l : Reenier. Amari. and Kersale . 



2002 ; Thompson. Kliem. and TorokL 201l[) . Close to the equator a bipolar struc 



ture emerges through the surface. The inner bulk has a positive current helicity. 
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Figure 14. Time series of a coronal ejection zoomed into the region of the ejection near the 
equator (6 = 7r/2), taken from Run A5. The dashed horizontal lines show the location of the 
surface at r = i?. Left column: normalized current helicity, /loRJ ■ B/{B^)t. Middle column: 
magnetic field, contours of r sin 6A^ are shown together with a color-scale representation of 
S^. The contours of r sinSA^ correspond to field lines of B in the r,9 plane, where solid lines 
represent clockwise magnetic field lines and the dashed ones counter-clockwise. Right column: 
density fluctuations Ap(t) = p{t) — {p)t- For all plots, the color-scale represents negative as 
dark blue and positive as light yellow. 
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Figure 15. X-point-like structure in the r,0 plane at the equator (0 = 7r/2) at t/T=2204 
zoomed into the ejection region, taken from Run A5. Contours of rAg are shown together 
with a color-scale representation of Bg; dark blue stands for negative and light yellow for 
positive values. The contours of rAg correspond to field lines of B in the r,<l> plane. The 
dashed horizontal lines show the location of the surface at r = i?. 



in Figure [T3] represented by a yellow color, and it pushes an arc with negative 
current helicity ahead of it; sec Figure 1141 Such bipolar ejections have been 
identified in earlier work (WBM) and compared with t he 'th ree-part structure' 
of coronal mass ejection, which is described in iLowl (|l996l ). The three parts 
consist of a prominence, which is similar to the bulk seen in our simulations, 
a front with an arc shaped structure corresponding to our arc, and a cavity 
between these two features. A bipola rity of tw i sted m agnetic field has also been 
seen in observed magnetic clouds bv lLi et~ai. (1982). Even though the domain 
of the simulation is larger in the direction than in WBM, the ejections arc 
much smaller, which is actually closer to the CMEs observed on the Sun. In 
the work of WBM the ejections have a size that corresponds to about 500 Mm, 
whereas in this work they seems to have a size corresponding to around 100 
Mm if scaled to the solar radius. The ejections seem to expand slightly, but 
no significant expansion rate can be measured using this resolution. Comparing 
with the forced turbulence runs, the difference in size is mostly due to the more 
complex and fluctuating magnetic field in convection runs. In the sequence of 
images of Figure 1141 an ejection near the equator reaches the outer boundary 
and leaves the domain. To investigate the mechanism driving the ejection, we 
look at the dynamics of the magnetic field in Figure where field lines of 
the azimuthally averaged mean field are shown as contours of rsin^^^, and 
colors represent together with the density fluctuations and current helicity. 
During the ejection, one may notice a strong concentration of magnetic field lines 
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Figure 16. Different properties of tlie ejection in tlie 6,0 plane at tlie surface (r = R) at 
i/r = 2204, taken from Run A5. Upper row, left panel: Contours of Ar are shown together with 
a color-scale representation of Br; dark blue stands for negative and light yellow for positive 
values. The contours of Ar correspond to field lines of that part of B that is solenoidal in 
the 6,rj) plane. Solid lines represent clockwise oriented magnetic field lines and dotted lines 
counter-clockwise ones. Middle panel: The arrows show {Ug,U^) and colors show Ur (blue 
corresponds to downflows). Right panel: Color-scale representation of the density p; dark blue 
stands for low and light yellow for high values. Lower row, left panel: Color-scale representation 
of specific entropy s; dark blue stands for low and light yellow for high values. Middle panel: 
Color-scale representation of the current density squared J^; dark blue stands for low and 
light yellow for high values. Right panel: current helicity J ■ B, color-scale as in Figure [131 The 
dashed line indicates the equator at 9 = it/2 = 90°. 



that are directed radially outwards. This concentration appears first beneath the 
surface and then emerges below the current helicity structure and follows it up 
into the coronal part. Investigating the direction of field lines of the mean field 
in the time series in Figure 1141 an X-point can be found. In the first panel, at 
r = 1.07 i? and 9 = 7r/2-|-0.1, the magnetic field lines form a junction-like shape. 
The dotted line represents a counter-clockwise oriented field loop, so at the two 
corners of the junction there are field lines with opposite signs. After around 14 
turnover times this "junction" has reconnected at the same position as where the 
ejection is detected. It appears that these two events are related to each other. 
Looking at the magnetic field line in the r,0 plane, which here is not averaged 
over the perpendicular direction (Figure USp , we identify a structure which has 
a shape similar to an X-point. 

The ejection causes also a strong variation in the density. If the time-averaged 
density profile is subtracted from instantaneous ones, the density fluctuations are 
obtained. After removing the density stratification one obtains Ap(t) = p{t) — 
{'p)t- We plot these density fluctuations, Ap{t), in the right column of Figure [TH 
to visualize the effect of the ejection on the density. The density in the ejection is 
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much lower than in the rest of the coronal part. However, the density variations 
are also associated with fluctuations in the specific entropy (As/cp « 0.01), 
which suggests that thermal buoyancy also plays a role. One interpretation could 
be that the strong magnetic field reduces the density to achieve total pressure 
equilibrium and the ejection rises partly because of magnetic buoyancy. Such an 
effect is also seen by inspecting other ejections. 

To characterize the emergence we plot different properties of the ejection in 
the 9,(1) plane; see Figure [TBI The magnetic field shows a strong concentration 
in its radial and azimutlial components. The concentration is associated with a 
downfiow in spite of it being a low-density region. It is interesting to note that in 
this case the gas velocity does not refiect the actual pattern speed. From the time 
evolution of the low-density region shown in Fig. I14[ we know that this region 
is moving radially upwards in a way that is consistent with a motion expected 
from buoyancy forces. In particular, the specific entropy has a high value in 
this region. In visualizations of the current density, we sec the formation of two 
current sheets. This leads to two current hclicity regions of opposite sign. 

When discussing coronal ejections, one is usually interested in the plasma 
/? parameter to characterize the corona. In our simplified coronal part, the 
plasma (3 does not decrease with radius, but it stays rather high, which is due 
to the low magnetic field strength, especially in the coronal part, even though 
B^j^s/Bpq =0.1-0.4 in the convection zone. The time-averaged value is always 
above 5 x 10**, and is therefore not comparable with the values in the solar 
corona, where the plasma /3 is very low because of the low density. There the 
magnetic field can drag dense plasma from the lower corona to its upper part. In 
our simulations the density stratification of the convection zone is much lower 
than in the Sun. Therefore, the density in the corona in our model is much higher 
and is closer to the density of the photosphere or the chromosphere. A rising 
magnetic fiux tube has formed a low-density region in its interior due to a higher 
magnetic pressure. As the tube rises further into the coronal part, the density 
inside the tube is still lower than that outside because the coronal density is 
rather high in our model. 

The simplification of a high plasma /? corona might not be suitable to de- 
scribe properly the mass fiux of the plasma drag ged by the magnetic field of 



the CME in the corona. Howe v er, th e early work of iMikic. Barnes, and Schnack 



(|l988l ). lOrtolani and Schnackl (|l993[ ). and IWiegelmannI (|2008l) has shown that 



an isothermal force-free approach (not to be confused with force-free magnetic 
equilibria) can describe the coronal magnetic field and even plasmoid ejections 
rather well. Note that in those papers the pressure gradient term was omitted, 
just like in the coronal part of WB. How important this really is remains unclear, 
because the pressure gradient term was not omitted in the work of WBM, which 
still showed ejections similar to those of WB. It would therefore be useful to 
compare our present model with one where the pressure gradient term is ignored 
in the coronal part, just like in WB. 

The ejection seen in Figures [TMTSl is not a single event — others follow in a 
recurrent fashion. However, the periodicity is not as clear as in previous work 
(WB.WBM). For Run A5, for example, we observe around five ejections during a 
time interval of about 1000 turnover times. A clearer indication for the recurrence 
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Figure 17. Recurrence of ejections shown by plotting the dependence of the dimcnsionless 
ratio /iqR J ■ B / {B^)t on time t/r and radius r in terms of the solar radius, taken from 
Run A5. The top panels show a narrow band in 8 in the northern hemisphere and the bottom 
ones in the southern hemisphere. We have also averaged in latitude from 4.1° to 19.5° (left 
panel) and 32.5° to 45.5° (right). Dark blue shades represent negative and light yellow positive 
values. The dashed horizontal lines show the location of the surface at r = R. Recurrence of 
ejections shown by plotting the dependence 



of the ejections can be seen in Figure [T71 where the normahzed current density 
is averaged over two narrow latitude bands in each hemisphere. The slope of 
structures in the outer parts in these rt diagrams gives an indication about the 
ejection speed Voj which turns out to be around one solar radius in 200-250 
turnover times. This translates to V^j/urms ~ 0.1, which is somewhat less than 
the values 0.2-0.5 found for the simulations of WBM. However, the mechanism 
which sets the time scale of ejections is at present still unclear. 

Given that gravity decreases with radius, there is i n principle the po ssibility 
of a radial wind with a critical point at = GM/2c^ ( Choudhuri . 19981 ). which 
would be at = 9.3 i?, i.e. well outside our coronal part. Because of this and 
the fact that we use closed boundary conditions with no mass flux out of the 
domain, no such wind can occur in our simulations. Using a boundary condition 
that would allow a mass flux in the radial direction could change the speed and 
the ejection properties significantly. Including a solar-like wind in a model can 
have two major effects, which require a much higher amount of computational 
resources. The radial variation of gravity applied in these simulations implies 
the presence of a critical point rather close to the surface of the convection zone. 
Therefore, if a wind were to develop, the resulting velocity in the convection 
zone would be too high for a dynamo to develop; the magnetic field would be 
blown out too quickly. Using instead a more realistic profile for the solar wind 
with a position of the critical point around = lOi?, the corresponding density 
stratification would be too strong to be stably resolved. 



4. Conclusions 

In the present paper we have presented an extension of the two-layer approach 
of WB and WBM by including a self-consistent rotating convection zone into the 
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model. Wc find a large-scale magnetic field generated by the convective turbulent 
motion in the convection zone. At moderate rotation rates, for a Coriolis number 
larger than 3, we obtain a differential rotation pattern showing super-rotation, 
i.e., an equator rotating faster than the poles. The dynamo solutions we find 
are different and some of them have a periodic oscillatory behavior, where the 
large-scale magnetic field does not change sign; only the strength is varying. At 
the maxima, the velocity is suppressed due to the backreaction via the Lorentz 
force. Small-scale magnetic structures seem to show an equatorward migration 
near the equator and a poleward one near the poles. 

Using a convectively driven dynamo complicates the generation of ejections 
into a coronal part due to lower relative kinetic hclicity. However, it was possible 
to produce ejections in two of the runs. The shape and the bipolar hclicity 
structure arc comparable with those of WBM. Due to the relatively high plasma 
/? in the outer parts of our model (compared with the solar corona), the ejections 
produce local minima of density which are carried along and ejected out to the 
top of the domain. The ejections occur recurrently, but not clearly periodically, 
which is similar to the Sun. 

Note that our results have to be interpreted cautiously, given the use of a 
simplistic solar atmosphere. We neglect the effects of high temperature and 
low plasma /?. However, we feel that the mechanism of emergence of magnetic 
structures driven by dynamo action from self-consistent convection may not 
strongly depend on these two conditions. This suggestion has to be proven in 
more detail in forthcoming work. 

An extension of the present work would require a detailed parameter study 
of cause and properties of the ejections. This also includes an advanced model 
for the solar corona with a lower plasma /3 and more efficient convection, which 
has a stronger stratification and is cooled by radiation. Another important as- 
pect would be the generation of a self-consistent solar wind which supports and 
interacts with the ejections. 
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